home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 7 / Apprentice-Release7.iso / Source Code / C / Applications / POV-Ray 3.0.2 / src / SOURCE / polysolv.c < prev    next >
Encoding:
C/C++ Source or Header  |  1997-05-18  |  29.8 KB  |  1,738 lines  |  [TEXT/MPS ]
  1. /****************************************************************************
  2. *                polysolv.c
  3. *
  4. *  This file was written by Alexander Enzmann.  He wrote the code for
  5. *  4th-6th order shapes and generously provided us these enhancements.
  6. *
  7. *  from Persistence of Vision(tm) Ray Tracer
  8. *  Copyright 1996 Persistence of Vision Team
  9. *---------------------------------------------------------------------------
  10. *  NOTICE: This source code file is provided so that users may experiment
  11. *  with enhancements to POV-Ray and to port the software to platforms other 
  12. *  than those supported by the POV-Ray Team.  There are strict rules under
  13. *  which you are permitted to use this file.  The rules are in the file
  14. *  named POVLEGAL.DOC which should be distributed with this file. If 
  15. *  POVLEGAL.DOC is not available or for more info please contact the POV-Ray
  16. *  Team Coordinator by leaving a message in CompuServe's Graphics Developer's
  17. *  Forum.  The latest version of POV-Ray may be found there as well.
  18. *
  19. * This program is based on the popular DKB raytracer version 2.12.
  20. * DKBTrace was originally written by David K. Buck.
  21. * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
  22. *
  23. *****************************************************************************/
  24.  
  25. #include "frame.h"
  26. #include "povray.h"
  27. #include "povproto.h"
  28. #include "vector.h"
  29. #include "polysolv.h"
  30.  
  31.  
  32.  
  33. /*****************************************************************************
  34. * Local preprocessor defines
  35. ******************************************************************************/
  36.  
  37. /*                  WARNING     WARNING    WARNING
  38.  
  39.    The following three constants have been defined so that quartic equations
  40.    will properly render on fpu/compiler combinations that only have 64 bit
  41.    IEEE precision.  Do not arbitrarily change any of these values.
  42.  
  43.    If you have a machine with a proper fpu/compiler combo - like a Mac with
  44.    a 68881, then use the native floating format (96 bits) and tune the
  45.    values for a little less tolerance: something like: factor1 = 1.0e15,
  46.    factor2 = -1.0e-7, factor3 = 1.0e-10.
  47.  
  48.    The meaning of the three constants are:
  49.       factor1 - the magnitude of difference between coefficients in a
  50.                 quartic equation at which the Sturmian root solver will
  51.                 kick in.  The Sturm solver is quite a bit slower than
  52.                 the algebraic solver, so the bigger this is, the faster
  53.                 the root solving will go.  The algebraic solver is less
  54.                 accurate so the bigger this is, the more likely you will
  55.                 get bad roots.
  56.  
  57.       factor2 - Tolerance value that defines how close the quartic equation
  58.                 is to being a square of a quadratic.  The closer this can
  59.                 get to zero before roots disappear, the less the chance
  60.                 you will get spurious roots.
  61.  
  62.       factor3 - Similar to factor2 at a later stage of the algebraic solver.
  63. */
  64.  
  65. #define FUDGE_FACTOR1 1.0e12
  66. #define FUDGE_FACTOR2 -1.0e-5
  67. #define FUDGE_FACTOR3 1.0e-7
  68.  
  69. /* Constants. */
  70.  
  71. #define TWO_M_PI_3  2.0943951023931954923084
  72. #define FOUR_M_PI_3 4.1887902047863909846168
  73.  
  74. /* Max number of iterations. */
  75.  
  76. #define MAX_ITERATIONS 50
  77.  
  78. /* A coefficient smaller than SMALL_ENOUGH is considered to be zero (0.0). */
  79.  
  80. #define SMALL_ENOUGH 1.0e-10
  81.  
  82. /* Smallest relative error we want. */
  83.  
  84. #define RELERROR 1.0e-12
  85.  
  86.  
  87.  
  88. /*****************************************************************************
  89. * Local typedefs
  90. ******************************************************************************/
  91.  
  92. typedef struct p
  93. {
  94.   int ord;
  95.   DBL coef[MAX_ORDER+1];
  96. }
  97. polynomial;
  98.  
  99.  
  100.  
  101. /*****************************************************************************
  102. * Local variables
  103. ******************************************************************************/
  104.  
  105.  
  106.  
  107. /*****************************************************************************
  108. * Static functions
  109. ******************************************************************************/
  110.  
  111. static int solve_quadratic PARAMS((DBL *x, DBL *y));
  112. static int solve_cubic PARAMS((DBL *x, DBL *y));
  113. static int solve_quartic PARAMS((DBL *x, DBL *y));
  114. static int polysolve PARAMS((int order, DBL *Coeffs, DBL *roots));
  115. static int modp PARAMS((polynomial *u, polynomial *v, polynomial *r));
  116. static int regula_falsa PARAMS((int order, DBL *coef, DBL a, DBL b, DBL *val));
  117. static int sbisect PARAMS((int np, polynomial *sseq, DBL min, DBL max, int atmin, int atmax, DBL *roots));
  118. static int numchanges PARAMS((int np, polynomial *sseq, DBL a));
  119. static DBL polyeval PARAMS((DBL x, int n, DBL *Coeffs));
  120. static int buildsturm PARAMS((int ord, polynomial *sseq));
  121. static int visible_roots PARAMS((int np, polynomial *sseq, int *atneg, int *atpos));
  122. static int difficult_coeffs PARAMS((int n, DBL *x));
  123.  
  124.  
  125.  
  126. /*****************************************************************************
  127. *
  128. * FUNCTION
  129. *
  130. *   modp
  131. *
  132. * INPUT
  133. *   
  134. * OUTPUT
  135. *   
  136. * RETURNS
  137. *   
  138. * AUTHOR
  139. *
  140. *   Alexander Enzmann
  141. *   
  142. * DESCRIPTION
  143. *
  144. *   Calculates the modulus of u(x) / v(x) leaving it in r.
  145. *   It returns 0 if r(x) is a constant.
  146. *
  147. *   NOTE: This function assumes the leading coefficient of v is 1 or -1.
  148. *
  149. * CHANGES
  150. *
  151. *   Okt 1996 : Added bug fix by Heiko Eissfeldt. [DB]
  152. *
  153. ******************************************************************************/
  154.  
  155. static int modp(u, v, r)
  156. polynomial *u, *v, *r;
  157. {
  158.   int k, j;
  159.  
  160.   *r = *u;
  161.  
  162.   if (v->coef[v->ord] < 0.0)
  163.   {
  164.     for (k = u->ord - v->ord - 1; k >= 0; k -= 2)
  165.     {
  166.       r->coef[k] = -r->coef[k];
  167.     }
  168.  
  169.     for (k = u->ord - v->ord; k >= 0; k--)
  170.     {
  171.       for (j = v->ord + k - 1; j >= k; j--)
  172.       {
  173.         r->coef[j] = -r->coef[j] - r->coef[v->ord + k] * v->coef[j - k];
  174.       }
  175.     }
  176.   }
  177.   else
  178.   {
  179.     for (k = u->ord - v->ord; k >= 0; k--)
  180.     {
  181.       for (j = v->ord + k - 1; j >= k; j--)
  182.       {
  183.         r->coef[j] -= r->coef[v->ord + k] * v->coef[j - k];
  184.       }
  185.     }
  186.   }
  187.  
  188.   k = v->ord - 1;
  189.  
  190.   while (k >= 0 && fabs(r->coef[k]) < SMALL_ENOUGH)
  191.   {
  192.     r->coef[k] = 0.0;
  193.  
  194.     k--;
  195.   }
  196.  
  197.   r->ord = (k < 0) ? 0 : k;
  198.  
  199.   return(r->ord);
  200. }
  201.  
  202.  
  203.  
  204. /*****************************************************************************
  205. *
  206. * FUNCTION
  207. *
  208. *   buildsturm
  209. *
  210. * INPUT
  211. *   
  212. * OUTPUT
  213. *   
  214. * RETURNS
  215. *   
  216. * AUTHOR
  217. *
  218. *   Alexander Enzmann
  219. *   
  220. * DESCRIPTION
  221. *
  222. *   Build the sturmian sequence for a polynomial.
  223. *
  224. * CHANGES
  225. *
  226. *   -
  227. *
  228. ******************************************************************************/
  229.  
  230. static int buildsturm(ord, sseq)
  231. int ord;
  232. polynomial *sseq;
  233. {
  234.   int i;
  235.   DBL f, *fp, *fc;
  236.   polynomial *sp;
  237.  
  238.   sseq[0].ord = ord;
  239.   sseq[1].ord = ord - 1;
  240.  
  241.   /* calculate the derivative and normalize the leading coefficient. */
  242.  
  243.   f = fabs(sseq[0].coef[ord] * ord);
  244.  
  245.   fp = sseq[1].coef;
  246.   fc = sseq[0].coef + 1;
  247.  
  248.   for (i = 1; i <= ord; i++)
  249.   {
  250.     *fp++ = *fc++ * i / f;
  251.   }
  252.  
  253.   /* construct the rest of the Sturm sequence */
  254.  
  255.   for (sp = sseq + 2; modp(sp - 2, sp - 1, sp); sp++)
  256.   {
  257.     /* reverse the sign and normalize */
  258.  
  259.     f = -fabs(sp->coef[sp->ord]);
  260.  
  261.     for (fp = &sp->coef[sp->ord]; fp >= sp->coef; fp--)
  262.     {
  263.       *fp /= f;
  264.     }
  265.   }
  266.  
  267.   /* reverse the sign */
  268.  
  269.   sp->coef[0] = -sp->coef[0];
  270.  
  271.   return(sp - sseq);
  272. }
  273.  
  274.  
  275.  
  276. /*****************************************************************************
  277. *
  278. * FUNCTION
  279. *
  280. *   visible_roots
  281. *
  282. * INPUT
  283. *   
  284. * OUTPUT
  285. *   
  286. * RETURNS
  287. *
  288. * AUTHOR
  289. *
  290. *   Alexander Enzmann
  291. *   
  292. * DESCRIPTION
  293. *
  294. *   Find out how many visible intersections there are.
  295. *
  296. * CHANGES
  297. *
  298. *   -
  299. *
  300. ******************************************************************************/
  301.  
  302. static int visible_roots(np, sseq, atzer, atpos)
  303. int np;
  304. polynomial *sseq;
  305. int *atzer, *atpos;
  306. {
  307.   int atposinf, atzero;
  308.   polynomial *s;
  309.   DBL f, lf;
  310.  
  311.   atposinf = atzero = 0;
  312.  
  313.   /* changes at positve infinity */
  314.  
  315.   lf = sseq[0].coef[sseq[0].ord];
  316.  
  317.   for (s = sseq + 1; s <= sseq + np; s++)
  318.   {
  319.     f = s->coef[s->ord];
  320.  
  321.     if (lf == 0.0 || lf * f < 0)
  322.     {
  323.       atposinf++;
  324.     }
  325.  
  326.     lf = f;
  327.   }
  328.  
  329.   /* Changes at zero */
  330.  
  331.   lf = sseq[0].coef[0];
  332.  
  333.   for (s = sseq+1; s <= sseq + np; s++)
  334.   {
  335.     f = s->coef[0];
  336.  
  337.     if (lf == 0.0 || lf * f < 0)
  338.     {
  339.       atzero++;
  340.     }
  341.  
  342.     lf = f;
  343.   }
  344.  
  345.   *atzer = atzero;
  346.   *atpos = atposinf;
  347.  
  348.   return(atzero - atposinf);
  349. }
  350.  
  351.  
  352.  
  353. /*****************************************************************************
  354. *
  355. * FUNCTION
  356. *
  357. *   numchanges
  358. *
  359. * INPUT
  360. *   
  361. * OUTPUT
  362. *   
  363. * RETURNS
  364. *   
  365. * AUTHOR
  366. *
  367. *   Alexander Enzmann
  368. *   
  369. * DESCRIPTION
  370. *
  371. *   Return the number of sign changes in the Sturm sequence in
  372. *   sseq at the value a.
  373. *
  374. * CHANGES
  375. *
  376. *   -
  377. *
  378. ******************************************************************************/
  379.  
  380. static int numchanges(np, sseq, a)
  381. int np;
  382. polynomial *sseq;
  383. DBL a;
  384. {
  385.   int changes;
  386.   DBL f, lf;
  387.   polynomial   *s;
  388.   changes = 0;
  389.  
  390.   lf = polyeval(a, sseq[0].ord, sseq[0].coef);
  391.  
  392.   for (s = sseq + 1; s <= sseq + np; s++)
  393.   {
  394.     f = polyeval(a, s->ord, s->coef);
  395.  
  396.     if (lf == 0.0 || lf * f < 0)
  397.     {
  398.       changes++;
  399.     }
  400.  
  401.     lf = f;
  402.   }
  403.  
  404.   return(changes);
  405. }
  406.  
  407.  
  408.  
  409. /*****************************************************************************
  410. *
  411. * FUNCTION
  412. *
  413. *   sbisect
  414. *
  415. * INPUT
  416. *   
  417. * OUTPUT
  418. *   
  419. * RETURNS
  420. *   
  421. * AUTHOR
  422. *
  423. *   Alexander Enzmann
  424. *   
  425. * DESCRIPTION
  426. *
  427. *   Uses a bisection based on the sturm sequence for the polynomial
  428. *   described in sseq to isolate intervals in which roots occur,
  429. *   the roots are returned in the roots array in order of magnitude.
  430. *
  431. *   NOTE: This routine has one severe bug: When the interval containing the
  432. *         root [min, max] has a root at one of its endpoints, as well as one
  433. *         within the interval, the root at the endpoint will be returned
  434. *         rather than the one inside.
  435. *
  436. * CHANGES
  437. *
  438. *   -
  439. *
  440. ******************************************************************************/
  441.  
  442. static int sbisect(np, sseq, min_value, max_value, atmin, atmax, roots)
  443. int np;
  444. polynomial *sseq;
  445. DBL min_value, max_value;
  446. int atmin, atmax;
  447. DBL *roots;
  448. {
  449.   DBL mid;
  450.   int n1, n2, its, atmid;
  451.   
  452.   if ((atmin - atmax) == 1)
  453.   {
  454.     /* first try using regula-falsa to find the root.  */
  455.  
  456.     if (regula_falsa(sseq->ord, sseq->coef, min_value, max_value, roots))
  457.     {
  458.       return(1);
  459.     }
  460.     else
  461.     {
  462.       /* That failed, so now find it by bisection */
  463.  
  464.       for (its = 0; its < MAX_ITERATIONS; its++)
  465.       {
  466.         mid = (min_value + max_value) / 2;
  467.  
  468.         atmid = numchanges(np, sseq, mid);
  469.         
  470.         /* The follow only happens if there is a bug.  And
  471.            unfortunately, there is. CEY 04/97 
  472.          */
  473.         if ((atmid<atmax) || (atmid>atmin))
  474.         {
  475.            return(0);
  476.         }
  477.  
  478.         if (fabs(mid) > RELERROR)
  479.         {
  480.           if (fabs((max_value - min_value) / mid) < RELERROR)
  481.           {
  482.             roots[0] = mid;
  483.  
  484.             return(1);
  485.           }
  486.         }
  487.         else
  488.         {
  489.           if (fabs(max_value - min_value) < RELERROR)
  490.           {
  491.             roots[0] = mid;
  492.  
  493.             return(1);
  494.           }
  495.         }
  496.  
  497.         if ((atmin - atmid) == 0)
  498.         {
  499.           min_value = mid;
  500.         }
  501.         else
  502.         {
  503.           max_value = mid;
  504.         }
  505.       }
  506.  
  507.       /* Bisection took too long - just return what we got */
  508.  
  509.       roots[0] = mid;
  510.  
  511.       return(1);
  512.     }
  513.   }
  514.  
  515.   /* There is more than one root in the interval.
  516.      Bisect to find new intervals. */
  517.  
  518.   for (its = 0; its < MAX_ITERATIONS; its++)
  519.   {
  520.     mid = (min_value + max_value) / 2;
  521.  
  522.     atmid = numchanges(np, sseq, mid);
  523.  
  524.     /* The follow only happens if there is a bug.  And
  525.        unfortunately, there is. CEY 04/97 
  526.      */
  527.     if ((atmid<atmax) || (atmid>atmin))
  528.     {
  529.        return(0);
  530.     }
  531.  
  532.     if (fabs(mid) > RELERROR)
  533.     {
  534.       if (fabs((max_value - min_value) / mid) < RELERROR)
  535.       {
  536.         roots[0] = mid;
  537.  
  538.         return(1);
  539.       }
  540.     }
  541.     else
  542.     {
  543.       if (fabs(max_value - min_value) < RELERROR)
  544.       {
  545.         roots[0] = mid;
  546.  
  547.         return(1);
  548.       }
  549.     }
  550.  
  551.     n1 = atmin - atmid;
  552.     n2 = atmid - atmax;
  553.  
  554.     if ((n1 != 0) && (n2 != 0))
  555.     {
  556.       n1 = sbisect(np, sseq, min_value, mid, atmin, atmid, roots);
  557.       n2 = sbisect(np, sseq, mid, max_value, atmid, atmax, &roots[n1]);
  558.  
  559.       return(n1 + n2);
  560.     }
  561.  
  562.     if (n1 == 0)
  563.     {
  564.       min_value = mid;
  565.     }
  566.     else
  567.     {
  568.       max_value = mid;
  569.     }
  570.   }
  571.  
  572.   /* Took too long to bisect.  Just return what we got. */
  573.  
  574.   roots[0] = mid;
  575.  
  576.   return(1);
  577. }
  578.  
  579.  
  580.  
  581. /*****************************************************************************
  582. *
  583. * FUNCTION
  584. *
  585. *   polyeval
  586. *
  587. * INPUT
  588. *   
  589. * OUTPUT
  590. *   
  591. * RETURNS
  592. *   
  593. * AUTHOR
  594. *
  595. *   Alexander Enzmann
  596. *   
  597. * DESCRIPTION
  598. *
  599. *   Evaluate the value of a polynomial at the given value x.
  600. *
  601. *   The coefficients are stored in c in the following order:
  602. *
  603. *     c[0] + c[1] * x + c[2] * x ^ 2 + c[3] * x ^ 3 + ...
  604. *
  605. * CHANGES
  606. *
  607. *   -
  608. *
  609. ******************************************************************************/
  610.  
  611. static DBL polyeval(x, n, Coeffs)
  612. DBL x, *Coeffs;
  613. int n;
  614. {
  615.   register int i;
  616.   DBL val;
  617.  
  618.   val = Coeffs[n];
  619.  
  620.   for (i = n-1; i >= 0; i--)
  621.   {
  622.     val = val * x + Coeffs[i];
  623.   }
  624.  
  625.   return(val);
  626. }
  627.  
  628.  
  629.  
  630. /*****************************************************************************
  631. *
  632. * FUNCTION
  633. *
  634. *   regular_falsa
  635. *
  636. * INPUT
  637. *   
  638. * OUTPUT
  639. *
  640. * RETURNS
  641. *
  642. * AUTHOR
  643. *
  644. *   Alexander Enzmann
  645. *
  646. * DESCRIPTION
  647. *
  648. *   Close in on a root by using regula-falsa.
  649. *
  650. * CHANGES
  651. *
  652. *   -
  653. *
  654. ******************************************************************************/
  655.  
  656. static int regula_falsa(order, coef, a, b, val)
  657. int order;
  658. DBL *coef;
  659. DBL a, b, *val;
  660. {
  661.   int its;
  662.   DBL fa, fb, x, fx, lfx;
  663.  
  664.   fa = polyeval(a, order, coef);
  665.   fb = polyeval(b, order, coef);
  666.  
  667.   if (fa * fb > 0.0)
  668.   {
  669.     return 0;
  670.   }
  671.  
  672.   if (fabs(fa) < SMALL_ENOUGH)
  673.   {
  674.     *val = a;
  675.  
  676.     return(1);
  677.   }
  678.  
  679.   if (fabs(fb) < SMALL_ENOUGH)
  680.   {
  681.     *val = b;
  682.  
  683.     return(1);
  684.   }
  685.  
  686.   lfx = fa;
  687.  
  688.   for (its = 0; its < MAX_ITERATIONS; its++)
  689.   {
  690.     x = (fb * a - fa * b) / (fb - fa);
  691.  
  692.     fx = polyeval(x, order, coef);
  693.  
  694.     if (fabs(x) > RELERROR)
  695.     {
  696.       if (fabs(fx / x) < RELERROR)
  697.       {
  698.         *val = x;
  699.  
  700.         return(1);
  701.       }
  702.     }
  703.     else
  704.     {
  705.       if (fabs(fx) < RELERROR)
  706.       {
  707.         *val = x;
  708.  
  709.         return(1);
  710.       }
  711.     }
  712.  
  713.     if (fa < 0)
  714.     {
  715.       if (fx < 0)
  716.       {
  717.         a = x;
  718.  
  719.         fa = fx;
  720.  
  721.         if ((lfx * fx) > 0)
  722.         {
  723.           fb /= 2;
  724.         }
  725.       }
  726.       else
  727.       {
  728.         b = x;
  729.  
  730.         fb = fx;
  731.  
  732.         if ((lfx * fx) > 0)
  733.         {
  734.           fa /= 2;
  735.         }
  736.       }
  737.     }
  738.     else
  739.     {
  740.       if (fx < 0)
  741.       {
  742.         b = x;
  743.  
  744.         fb = fx;
  745.  
  746.         if ((lfx * fx) > 0)
  747.         {
  748.           fa /= 2;
  749.         }
  750.       }
  751.       else
  752.       {
  753.         a = x;
  754.  
  755.         fa = fx;
  756.  
  757.         if ((lfx * fx) > 0)
  758.         {
  759.           fb /= 2;
  760.         }
  761.       }
  762.     }
  763.  
  764.     /* Check for underflow in the domain */
  765.  
  766.     if (fabs(b-a) < RELERROR)
  767.     {
  768.       *val = x;
  769.  
  770.       return(1);
  771.     }
  772.  
  773.     lfx = fx;
  774.   }
  775.  
  776.   return(0);
  777. }
  778.  
  779.  
  780.  
  781. /*****************************************************************************
  782. *
  783. * FUNCTION
  784. *
  785. *   solve_quadratic
  786. *
  787. * INPUT
  788. *
  789. * OUTPUT
  790. *
  791. * RETURNS
  792. *   
  793. * AUTHOR
  794. *
  795. *   Alexander Enzmann
  796. *   
  797. * DESCRIPTION
  798. *
  799. *   Solve the quadratic equation:
  800. *
  801. *     x[0] * x^2 + x[1] * x + x[2] = 0.
  802. *
  803. *   The value returned by this function is the number of real roots.
  804. *   The roots themselves are returned in y[0], y[1].
  805. *
  806. * CHANGES
  807. *
  808. *   -
  809. *
  810. ******************************************************************************/
  811.  
  812. static int solve_quadratic(x, y)
  813. DBL *x, *y;
  814. {
  815.   DBL d, t, a, b, c;
  816.  
  817.   a = x[0];
  818.   b = -x[1];
  819.   c = x[2];
  820.  
  821.   if (a == 0.0)
  822.   {
  823.     if (b == 0.0)
  824.     {
  825.       return(0);
  826.     }
  827.  
  828.     y[0] = c / b;
  829.  
  830.     return(1);
  831.   }
  832.  
  833.   d = b * b - 4.0 * a * c;
  834.  
  835.   /* Treat values of d around 0 as 0. */
  836.  
  837.   if ((d > -SMALL_ENOUGH) && (d < SMALL_ENOUGH))
  838.   {
  839.     y[0] = 0.5 * b / a;
  840.  
  841.     return(1);
  842.   }
  843.   else
  844.   {
  845.     if (d < 0.0)
  846.     {
  847.       return(0);
  848.     }
  849.   }
  850.  
  851.   d = sqrt(d);
  852.  
  853.   t = 2.0 * a;
  854.  
  855.   y[0] = (b + d) / t;
  856.   y[1] = (b - d) / t;
  857.  
  858.   return(2);
  859. }
  860.  
  861.  
  862.  
  863. /*****************************************************************************
  864. *
  865. * FUNCTION
  866. *
  867. *   solve_cubic
  868. *
  869. * INPUT
  870. *   
  871. * OUTPUT
  872. *   
  873. * RETURNS
  874. *   
  875. * AUTHOR
  876. *
  877. *   Alexander Enzmann
  878. *   
  879. * DESCRIPTION
  880. *
  881. *
  882. *   Solve the cubic equation:
  883. *
  884. *     x[0] * x^3 + x[1] * x^2 + x[2] * x + x[3] = 0.
  885. *
  886. *   The result of this function is an integer that tells how many real
  887. *   roots exist.  Determination of how many are distinct is up to the
  888. *   process that calls this routine.  The roots that exist are stored
  889. *   in (y[0], y[1], y[2]).
  890. *
  891. *   NOTE: This function relies very heavily on trigonometric functions and
  892. *         the square root function.  If an alternative solution is found
  893. *         that does not rely on transcendentals this code will be replaced.
  894. *
  895. * CHANGES
  896. *
  897. *   -
  898. *
  899. ******************************************************************************/
  900.  
  901. static int solve_cubic(x, y)
  902. DBL *x, *y;
  903. {
  904.   DBL Q, R, Q3, R2, sQ, d, an, theta;
  905.   DBL A2, a0, a1, a2, a3;
  906.  
  907.   a0 = x[0];
  908.  
  909.   if (a0 == 0.0)
  910.   {
  911.     return(solve_quadratic(&x[1], y));
  912.   }
  913.   else
  914.   {
  915.     if (a0 != 1.0)
  916.     {
  917.       a1 = x[1] / a0;
  918.       a2 = x[2] / a0;
  919.       a3 = x[3] / a0;
  920.     }
  921.     else
  922.     {
  923.       a1 = x[1];
  924.       a2 = x[2];
  925.       a3 = x[3];
  926.     }
  927.   }
  928.  
  929.   A2 = a1 * a1;
  930.  
  931.   Q = (A2 - 3.0 * a2) / 9.0;
  932.  
  933.   /* Modified to save some multiplications and to avoid a floating point
  934.      exception that occured with DJGPP and full optimization. [DB 8/94] */
  935.  
  936.   R = (a1 * (A2 - 4.5 * a2) + 13.5 * a3) / 27.0;
  937.  
  938.   Q3 = Q * Q * Q;
  939.  
  940.   R2 = R * R;
  941.  
  942.   d = Q3 - R2;
  943.  
  944.   an = a1 / 3.0;
  945.  
  946.   if (d >= 0.0)
  947.   {
  948.     /* Three real roots. */
  949.  
  950.     d = R / sqrt(Q3);
  951.  
  952.     theta = acos(d) / 3.0;
  953.  
  954.     sQ = -2.0 * sqrt(Q);
  955.  
  956.     y[0] = sQ * cos(theta) - an;
  957.     y[1] = sQ * cos(theta + TWO_M_PI_3) - an;
  958.     y[2] = sQ * cos(theta + FOUR_M_PI_3) - an;
  959.  
  960.     return(3);
  961.   }
  962.   else
  963.   {
  964.     sQ = pow(sqrt(R2 - Q3) + fabs(R), 1.0 / 3.0);
  965.  
  966.     if (R < 0)
  967.     {
  968.       y[0] = (sQ + Q / sQ) - an;
  969.     }
  970.     else
  971.     {
  972.       y[0] = -(sQ + Q / sQ) - an;
  973.     }
  974.  
  975.     return(1);
  976.   }
  977. }
  978.  
  979. #ifdef USE_NEW_DIFFICULT_COEFFS
  980.  
  981. /*****************************************************************************
  982. *
  983. * FUNCTION
  984. *
  985. *   difficult_coeffs
  986. *
  987. * INPUT
  988. *   
  989. * OUTPUT
  990. *   
  991. * RETURNS
  992. *   
  993. * AUTHOR
  994. *
  995. *   Alexander Enzmann
  996. *
  997. * DESCRIPTION
  998. *
  999. *   Test to see if any coeffs are more than 6 orders of magnitude
  1000. *   larger than the smallest.
  1001. *
  1002. * CHANGES
  1003. *
  1004. *   -
  1005. *
  1006. ******************************************************************************/
  1007.  
  1008. static int difficult_coeffs(n, x)
  1009. int n;
  1010. DBL *x;
  1011. {
  1012.   int i, flag = 0;
  1013.   DBL t, biggest;
  1014.  
  1015.   biggest = fabs(x[0]);
  1016.  
  1017.   for (i = 1; i <= n; i++)
  1018.   {
  1019.     t = fabs(x[i]);
  1020.  
  1021.     if (t > biggest)
  1022.     {
  1023.       biggest = t;
  1024.     }
  1025.   }
  1026.  
  1027.   /* Everything is zero no sense in doing any more */
  1028.  
  1029.   if (biggest == 0.0)
  1030.   {
  1031.     return(0);
  1032.   }
  1033.  
  1034.   for (i = 0; i <= n; i++)
  1035.   {
  1036.     if (x[i] != 0.0)
  1037.     {
  1038.       if (fabs(biggest / x[i]) > FUDGE_FACTOR1)
  1039.       {
  1040.         x[i] = 0.0;
  1041.         flag = 1;
  1042.       }
  1043.     }
  1044.   }
  1045.  
  1046.   return(flag);
  1047. }
  1048. #else
  1049. /*****************************************************************************
  1050. *
  1051. * FUNCTION
  1052. *
  1053. *   difficult_coeffs
  1054. *
  1055. * INPUT
  1056. *   
  1057. * OUTPUT
  1058. *   
  1059. * RETURNS
  1060. *   
  1061. * AUTHOR
  1062. *
  1063. *   Alexander Enzmann
  1064. *   
  1065. * DESCRIPTION
  1066. *
  1067. *   Test to see if any coeffs are more than 6 orders of magnitude
  1068. *   larger than the smallest (function from POV 1.0) [DB 8/94].
  1069. *
  1070. * CHANGES
  1071. *
  1072. *   -
  1073. *
  1074. ******************************************************************************/
  1075.  
  1076. static int difficult_coeffs(n, x)
  1077. int n;
  1078. DBL *x;
  1079. {
  1080.   int i;
  1081.   DBL biggest;
  1082.  
  1083.   biggest = 0.0;
  1084.  
  1085.   for (i = 0; i <= n; i++)
  1086.   {
  1087.     if (fabs(x[i]) > biggest)
  1088.     {
  1089.       biggest = x[i];
  1090.     }
  1091.   }
  1092.  
  1093.   /* Everything is zero no sense in doing any more */
  1094.  
  1095.   if (biggest == 0.0)
  1096.   {
  1097.     return(0);
  1098.   }
  1099.  
  1100.   for (i = 0; i <= n; i++)
  1101.   {
  1102.     if (x[i] != 0.0)
  1103.     {
  1104.       if (fabs(biggest / x[i]) > FUDGE_FACTOR1)
  1105.       {
  1106.         return(1);
  1107.       }
  1108.     }
  1109.   }
  1110.  
  1111.   return(0);
  1112. }
  1113.  
  1114. #endif
  1115.  
  1116. #ifdef TEST_SOLVER
  1117. /*****************************************************************************
  1118. *
  1119. * FUNCTION
  1120. *
  1121. *   solve_quartic
  1122. *
  1123. * INPUT
  1124. *   
  1125. * OUTPUT
  1126. *   
  1127. * RETURNS
  1128. *   
  1129. * AUTHOR
  1130. *
  1131. *   Alexander Enzmann
  1132. *
  1133. * DESCRIPTION
  1134. *
  1135. *   The old way of solving quartics algebraically.
  1136. *   This is an adaptation of the method of Lodovico Ferrari (Circa 1731).
  1137. *
  1138. * CHANGES
  1139. *
  1140. *   -
  1141. *
  1142. ******************************************************************************/
  1143.  
  1144. static int solve_quartic(x, results)
  1145. DBL *x, *results;
  1146. {
  1147.   DBL cubic[4], roots[3];
  1148.   DBL a0, a1, y, d1, x1, t1, t2;
  1149.   DBL c0, c1, c2, c3, c4, d2, q1, q2;
  1150.   int i;
  1151.  
  1152.   c0 = x[0];
  1153.  
  1154.   if (c0 != 1.0)
  1155.   {
  1156.     c1 = x[1] / c0;
  1157.     c2 = x[2] / c0;
  1158.     c3 = x[3] / c0;
  1159.     c4 = x[4] / c0;
  1160.   }
  1161.   else
  1162.   {
  1163.     c1 = x[1];
  1164.     c2 = x[2];
  1165.     c3 = x[3];
  1166.     c4 = x[4];
  1167.   }
  1168.  
  1169.   /* The first step is to take the original equation:
  1170.  
  1171.        x^4 + b*x^3 + c*x^2 + d*x + e = 0
  1172.  
  1173.      and rewrite it as:
  1174.  
  1175.        x^4 + b*x^3 = -c*x^2 - d*x - e,
  1176.  
  1177.      adding (b*x/2)^2 + (x^2 + b*x/2)y + y^2/4 to each side gives a
  1178.      perfect square on the lhs:
  1179.  
  1180.        (x^2 + b*x/2 + y/2)^2 = (b^2/4 - c + y)x^2 + (b*y/2 - d)x + y^2/4 - e
  1181.  
  1182.      By choosing the appropriate value for y, the rhs can be made a perfect
  1183.      square also.  This value is found when the rhs is treated as a quadratic
  1184.      in x with the discriminant equal to 0.  This will be true when:
  1185.  
  1186.        (b*y/2 - d)^2 - 4.0 * (b^2/4 - c*y)*(y^2/4 - e) = 0, or
  1187.  
  1188.        y^3 - c*y^2 + (b*d - 4*e)*y - b^2*e + 4*c*e - d^2 = 0.
  1189.  
  1190.      This is called the resolvent of the quartic equation.  */
  1191.  
  1192.   a0 = 4.0 * c4;
  1193.  
  1194.   cubic[0] = 1.0;
  1195.   cubic[1] = -1.0 * c2;
  1196.   cubic[2] = c1 * c3 - a0;
  1197.   cubic[3] = a0 * c2 - c1 * c1 * c4 - c3 * c3;
  1198.  
  1199.   i = solve_cubic(&cubic[0], &roots[0]);
  1200.  
  1201.   if (i > 0)
  1202.   {
  1203.     y = roots[0];
  1204.   }
  1205.   else
  1206.   {
  1207.     return(0);
  1208.   }
  1209.  
  1210.   /* What we are left with is a quadratic squared on the lhs and a
  1211.      linear term on the right.  The linear term has one of two signs,
  1212.      take each and add it to the lhs.  The form of the quartic is now:
  1213.  
  1214.        a' = b^2/4 - c + y,    b' = b*y/2 - d, (from rhs quadritic above)
  1215.  
  1216.        (x^2 + b*x/2 + y/2) = +sqrt(a'*(x + 1/2 * b'/a')^2), and
  1217.        (x^2 + b*x/2 + y/2) = -sqrt(a'*(x + 1/2 * b'/a')^2).
  1218.  
  1219.      By taking the linear term from each of the right hand sides and
  1220.      adding to the appropriate part of the left hand side, two quadratic
  1221.      formulas are created.  By solving each of these the four roots of
  1222.      the quartic are determined.
  1223.   */
  1224.  
  1225.   i = 0;
  1226.  
  1227.   a0 = c1 / 2.0;
  1228.   a1 = y / 2.0;
  1229.  
  1230.   t1 = a0 * a0 - c2 + y;
  1231.  
  1232.   if (t1 < 0.0)
  1233.   {
  1234.     if (t1 > FUDGE_FACTOR2)
  1235.     {
  1236.       t1 = 0.0;
  1237.     }
  1238.     else
  1239.     {
  1240.       /* First Special case, a' < 0 means all roots are complex. */
  1241.  
  1242.       return(0);
  1243.       }
  1244.     }
  1245.   }
  1246.  
  1247.   if (t1 < FUDGE_FACTOR3)
  1248.   {
  1249.     /* Second special case, the "x" term on the right hand side above
  1250.        has vanished.  In this case:
  1251.  
  1252.          (x^2 + b*x/2 + y/2) = +sqrt(y^2/4 - e), and
  1253.          (x^2 + b*x/2 + y/2) = -sqrt(y^2/4 - e).  */
  1254.  
  1255.     t2 = a1 * a1 - c4;
  1256.  
  1257.     if (t2 < 0.0)
  1258.     {
  1259.       return(0);
  1260.     }
  1261.  
  1262.     x1 = 0.0;
  1263.     d1 = sqrt(t2);
  1264.   }
  1265.   else
  1266.   {
  1267.     x1 = sqrt(t1);
  1268.     d1 = 0.5 * (a0 * y - c3) / x1;
  1269.   }
  1270.  
  1271.   /* Solve the first quadratic */
  1272.  
  1273.   q1 = -a0 - x1;
  1274.   q2 = a1 + d1;
  1275.   d2 = q1 * q1 - 4.0 * q2;
  1276.  
  1277.   if (d2 >= 0.0)
  1278.   {
  1279.     d2 = sqrt(d2);
  1280.  
  1281.     results[0] = 0.5 * (q1 + d2);
  1282.     results[1] = 0.5 * (q1 - d2);
  1283.  
  1284.     i = 2;
  1285.   }
  1286.  
  1287.   /* Solve the second quadratic */
  1288.  
  1289.   q1 = q1 + x1 + x1;
  1290.   q2 = a1 - d1;
  1291.   d2 = q1 * q1 - 4.0 * q2;
  1292.  
  1293.   if (d2 >= 0.0)
  1294.   {
  1295.     d2 = sqrt(d2);
  1296.  
  1297.     results[i++] = 0.5 * (q1 + d2);
  1298.     results[i++] = 0.5 * (q1 - d2);
  1299.   }
  1300.  
  1301.   return(i);
  1302. }
  1303. #else
  1304. /*****************************************************************************
  1305. *
  1306. * FUNCTION
  1307. *
  1308. *   solve_quartic
  1309. *
  1310. * INPUT
  1311. *   
  1312. * OUTPUT
  1313. *   
  1314. * RETURNS
  1315. *   
  1316. * AUTHOR
  1317. *
  1318. *   Alexander Enzmann
  1319. *
  1320. * DESCRIPTION
  1321. *
  1322. *   Solve a quartic using the method of Francois Vieta (Circa 1735).
  1323. *
  1324. * CHANGES
  1325. *
  1326. *   -
  1327. *
  1328. ******************************************************************************/
  1329.  
  1330. static int solve_quartic(x, results)
  1331. DBL *x, *results;
  1332. {
  1333.   DBL cubic[4], roots[3];
  1334.   DBL c12, z, p, q, q1, q2, r, d1, d2;
  1335.   DBL c0, c1, c2, c3, c4;
  1336.   int i;
  1337.  
  1338.   /* Make sure the quartic has a leading coefficient of 1.0 */
  1339.  
  1340.   c0 = x[0];
  1341.  
  1342.   if (c0 != 1.0)
  1343.   {
  1344.     c1 = x[1] / c0;
  1345.     c2 = x[2] / c0;
  1346.     c3 = x[3] / c0;
  1347.     c4 = x[4] / c0;
  1348.   }
  1349.   else
  1350.   {
  1351.     c1 = x[1];
  1352.     c2 = x[2];
  1353.     c3 = x[3];
  1354.     c4 = x[4];
  1355.   }
  1356.  
  1357.   /* Compute the cubic resolvant */
  1358.  
  1359.   c12 = c1 * c1;
  1360.   p = -0.375 * c12 + c2;
  1361.   q = 0.125 * c12 * c1 - 0.5 * c1 * c2 + c3;
  1362.   r = -0.01171875 * c12 * c12 + 0.0625 * c12 * c2 - 0.25 * c1 * c3 + c4;
  1363.  
  1364.   cubic[0] = 1.0;
  1365.   cubic[1] = -0.5 * p;
  1366.   cubic[2] = -r;
  1367.   cubic[3] = 0.5 * r * p - 0.125 * q * q;
  1368.  
  1369.   i = solve_cubic(cubic, roots);
  1370.  
  1371.   if (i > 0)
  1372.   {
  1373.     z = roots[0];
  1374.   }
  1375.   else
  1376.   {
  1377.     return(0);
  1378.   }
  1379.  
  1380.   d1 = 2.0 * z - p;
  1381.  
  1382.   if (d1 < 0.0)
  1383.   {
  1384.     if (d1 > -SMALL_ENOUGH)
  1385.     {
  1386.       d1 = 0.0;
  1387.     }
  1388.     else
  1389.     {
  1390.       return(0);
  1391.     }
  1392.   }
  1393.  
  1394.   if (d1 < SMALL_ENOUGH)
  1395.   {
  1396.     d2 = z * z - r;
  1397.  
  1398.     if (d2 < 0.0)
  1399.     {
  1400.       return(0);
  1401.     }
  1402.  
  1403.     d2 = sqrt(d2);
  1404.   }
  1405.   else
  1406.   {
  1407.     d1 = sqrt(d1);
  1408.     d2 = 0.5 * q / d1;
  1409.   }
  1410.  
  1411.   /* Set up useful values for the quadratic factors */
  1412.  
  1413.   q1 = d1 * d1;
  1414.   q2 = -0.25 * c1;
  1415.  
  1416.   i = 0;
  1417.  
  1418.   /* Solve the first quadratic */
  1419.  
  1420.   p = q1 - 4.0 * (z - d2);
  1421.  
  1422.   if (p == 0)
  1423.   {
  1424.     results[i++] = -0.5 * d1 - q2;
  1425.   }
  1426.   else
  1427.   {
  1428.     if (p > 0)
  1429.     {
  1430.       p = sqrt(p);
  1431.       results[i++] = -0.5 * (d1 + p) + q2;
  1432.       results[i++] = -0.5 * (d1 - p) + q2;
  1433.     }
  1434.   }
  1435.  
  1436.   /* Solve the second quadratic */
  1437.  
  1438.   p = q1 - 4.0 * (z + d2);
  1439.  
  1440.   if (p == 0)
  1441.   {
  1442.     results[i++] = 0.5 * d1 - q2;
  1443.   }
  1444.   else
  1445.   {
  1446.     if (p > 0)
  1447.     {
  1448.       p = sqrt(p);
  1449.       results[i++] = 0.5 * (d1 + p) + q2;
  1450.       results[i++] = 0.5 * (d1 - p) + q2;
  1451.     }
  1452.   }
  1453.  
  1454.   return(i);
  1455. }
  1456. #endif
  1457.  
  1458.  
  1459.  
  1460. /*****************************************************************************
  1461. *
  1462. * FUNCTION
  1463. *
  1464. *   polysolve
  1465. *
  1466. * INPUT
  1467. *   
  1468. * OUTPUT
  1469. *   
  1470. * RETURNS
  1471. *   
  1472. * AUTHOR
  1473. *
  1474. *   Alexander Enzmann
  1475. *   
  1476. * DESCRIPTION
  1477. *
  1478. *   Root solver based on the Sturm sequences for a polynomial.
  1479. *
  1480. * CHANGES
  1481. *
  1482. *   Okt 1996 : Added bug fix by Heiko Eissfeldt. [DB]
  1483. *
  1484. ******************************************************************************/
  1485.  
  1486. static int polysolve(order, Coeffs, roots)
  1487. int order;
  1488. DBL *Coeffs, *roots;
  1489. {
  1490.   polynomial sseq[MAX_ORDER+1];
  1491.   DBL min_value, max_value;
  1492.   int i, nroots, np, atmin, atmax;
  1493.  
  1494.   /* Put the coefficients into the top of the stack. */
  1495.  
  1496.   for (i = 0; i <= order; i++)
  1497.   {
  1498.     sseq[0].coef[order-i] = Coeffs[i] / Coeffs[0];
  1499.   }
  1500.  
  1501.   /* Build the Sturm sequence */
  1502.  
  1503.   np = buildsturm(order, &sseq[0]);
  1504.  
  1505.  
  1506.   /* Get the total number of visible roots */
  1507.  
  1508.   if ((nroots = visible_roots(np, sseq, &atmin, &atmax)) == 0)
  1509.   {
  1510.     return(0);
  1511.   }
  1512.  
  1513.   /* Bracket the roots */
  1514.  
  1515.   min_value = 0.0;
  1516.   max_value = Max_Distance;
  1517.  
  1518.   atmin = numchanges(np, sseq, min_value);
  1519.   atmax = numchanges(np, sseq, max_value);
  1520.  
  1521.   nroots = atmin - atmax;
  1522.  
  1523.   if (nroots == 0)
  1524.   {
  1525.     return(0);
  1526.   }
  1527.  
  1528.   /* perform the bisection. */
  1529.  
  1530.   return(sbisect(np, sseq, min_value, max_value, atmin, atmax, roots));
  1531. }
  1532.  
  1533.  
  1534.  
  1535. /*****************************************************************************
  1536. *
  1537. * FUNCTION
  1538. *
  1539. *   Solve_Polynomial
  1540. *
  1541. * INPUT
  1542. *
  1543. *   n       - order of polynomial
  1544. *   c       - coefficients
  1545. *   r       - roots
  1546. *   sturm   - TRUE, if sturm should be used for n=3,4
  1547. *   epsilon - Tolerance to discard small root
  1548. *   
  1549. * OUTPUT
  1550. *
  1551. *   r
  1552. *   
  1553. * RETURNS
  1554. *
  1555. *   int - number of roots found
  1556. *   
  1557. * AUTHOR
  1558. *
  1559. *   Dieter Bayer
  1560. *   
  1561. * DESCRIPTION
  1562. *
  1563. *   Solve the polynomial equation
  1564. *
  1565. *     c[0] * x ^ n + c[1] * x ^ (n-1) + ... + c[n-1] * x + c[n] = 0
  1566. *
  1567. *   If the equation has a root r, |r| < epsilon, this root is eliminated
  1568. *   and the equation of order n-1 will be solved. This will avoid the problem
  1569. *   of "surface acne" in (most) cases while increasing the speed of the
  1570. *   root solving (polynomial's order is reduced by one).
  1571. *
  1572. *   WARNING: This function can only be used for polynomials if small roots
  1573. *   (i.e. |x| < epsilon) are not needed. This is the case for ray/object
  1574. *   intersection tests because only intersecions with t > 0 are valid.
  1575. *
  1576. *   NOTE: Only one root at x = 0 will be eliminated.
  1577. *
  1578. *   NOTE: If epsilon = 0 no roots will be eliminated.
  1579. *
  1580. *
  1581. *   The method and idea for root elimination was taken from:
  1582. *
  1583. *     Han-Wen Nienhuys, "Polynomials", Ray Tracing News, July 6, 1994,
  1584. *     Volume 7, Number 3
  1585. *
  1586. *
  1587. * CHANGES
  1588. *
  1589. *   Jul 1994 : Creation.
  1590. *
  1591. ******************************************************************************/
  1592.  
  1593. int Solve_Polynomial(n, c0, r, sturm, epsilon)
  1594. int n, sturm;
  1595. DBL *c0, *r, epsilon;
  1596. {
  1597.   int roots, i;
  1598.   DBL *c;
  1599.  
  1600.   Increase_Counter(stats[Polynomials_Tested]);
  1601.  
  1602.   roots = 0;
  1603.  
  1604.   /*
  1605.    * Determine the "real" order of the polynomial, i.e.
  1606.    * eliminate small leading coefficients.
  1607.    */
  1608.  
  1609.   i = 0;
  1610.  
  1611.   while ((fabs(c0[i]) < SMALL_ENOUGH) && (i < n))
  1612.   {
  1613.     i++;
  1614.   }
  1615.  
  1616.   n -= i;
  1617.  
  1618.   c = &c0[i];
  1619.  
  1620.   switch (n)
  1621.   {
  1622.     case 0:
  1623.  
  1624.       break;
  1625.  
  1626.     case 1:
  1627.  
  1628.       /* Solve linear polynomial. */
  1629.  
  1630.       if (c[0] != 0.0)
  1631.       {
  1632.         r[roots++] = -c[1] / c[0];
  1633.       }
  1634.  
  1635.       break;
  1636.  
  1637.     case 2:
  1638.  
  1639.       /* Solve quadratic polynomial. */
  1640.  
  1641.       roots = solve_quadratic(c, r);
  1642.  
  1643.       break;
  1644.  
  1645.     case 3:
  1646.  
  1647.       /* Root elimination? */
  1648.  
  1649.       if (epsilon > 0.0)
  1650.       {
  1651.         if ((c[2] != 0.0) && (fabs(c[3]/c[2]) < epsilon))
  1652.         {
  1653.           Increase_Counter(stats[Roots_Eliminated]);
  1654.  
  1655.           roots = solve_quadratic(c, r);
  1656.  
  1657.           break;
  1658.         }
  1659.       }
  1660.  
  1661.       /* Solve cubic polynomial. */
  1662.  
  1663.       if (sturm)
  1664.       {
  1665.         roots = polysolve(3, c, r);
  1666.       }
  1667.       else
  1668.       {
  1669.         roots = solve_cubic(c, r);
  1670.       }
  1671.  
  1672.       break;
  1673.  
  1674.     case 4:
  1675.  
  1676.       /* Root elimination? */
  1677.  
  1678.       if (epsilon > 0.0)
  1679.       {
  1680.         if ((c[3] != 0.0) && (fabs(c[4]/c[3]) < epsilon))
  1681.         {
  1682.           Increase_Counter(stats[Roots_Eliminated]);
  1683.  
  1684.           if (sturm)
  1685.           {
  1686.             roots = polysolve(3, c, r);
  1687.           }
  1688.           else
  1689.           {
  1690.             roots = solve_cubic(c, r);
  1691.           }
  1692.  
  1693.           break;
  1694.         }
  1695.       }
  1696.  
  1697.       /* Test for difficult coeffs. */
  1698.  
  1699.       if (difficult_coeffs(4, c))
  1700.       {
  1701.         sturm = TRUE;
  1702.       }
  1703.  
  1704.       /* Solve quartic polynomial. */
  1705.  
  1706.       if (sturm)
  1707.       {
  1708.         roots = polysolve(4, c, r);
  1709.       }
  1710.       else
  1711.       {
  1712.         roots = solve_quartic(c, r);
  1713.       }
  1714.  
  1715.       break;
  1716.  
  1717.     default:
  1718.  
  1719.       if (epsilon > 0.0)
  1720.       {
  1721.         if ((c[n-1] != 0.0) && (fabs(c[n]/c[n-1]) < epsilon))
  1722.         {
  1723.           Increase_Counter(stats[Roots_Eliminated]);
  1724.  
  1725.           roots = polysolve(n-1, c, r);
  1726.         }
  1727.       }
  1728.  
  1729.       /* Solve n-th order polynomial. */
  1730.  
  1731.       roots = polysolve(n, c, r);
  1732.  
  1733.       break;
  1734.   }
  1735.  
  1736.   return(roots);
  1737. }
  1738.